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Abstract 

Previous  one-dimensional  Hall  thruster  models 
predict  some  observed  features  correctly,  but  do  not 
adequately  treat  the  losses  of  energy  and  ionization 
to  the  side  walls  of  the  channel.  Because  those 
losses  have  important  effects  on  the  discharge 
dynamics  and  the  thruster’s  efficiency  and  lifetime, 
their  omission  severely  limits  the  utility  of  a  model. 
The  model  reported  here  is  obtained  from  a  two- 
dimensional  description  by  solving  analytically  for 
the  radial  variations  of  densities  and  velocities  and 
then  averaging  over  the  radial  coordinate  in  a  way 
that  retains  the  effects  of  the  side  walls.  Our  model 
also  includes  ion  pressure,  second  ionization,  a 
varying  electron  temperature,  a  diverging  channel 
and  a  non-radial  B  field.  It  calculates  the  ion  flux 
into  the  walls  as  well  as  the  incident  ion  angles  and 
energies  as  functions  of  axial  position  to  predict  the 
profile  of  erosion  of  the  walls  by  sputtering.  The 
collision  and  ionization  rates  are  all  treated  as 
functions  of  the  electron  temperature.  The 
derivation  of  the  equations  is  explained  and  typical 
results  of  numerical  solutions  are  presented. 

1.0  Introduction 

Hall  thrusters  have  been  modeled  numerically  using 
1-D  fluid  equations11'63  and  2-D  hybrid  and  particle 
descriptions.17'141  To  include  enough  particles  for 
accurate  predictions,  the  2-D  codes  require  long 
computations.  The  simpler  1-D  models  run  faster 
and  can  therefore  survey  more  cases,  but  the  utility 
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of  existing  1-D  models  is  limited  by  their  lack  of  an 
accurate  treatment  of  the  effects  of  the  side  walls  of 
the  channel,  which  are  a  sink  for  ionization  and 
energy.  Plasma  recombines  at  the  walls  and  the 
discharge  electrons  are  cooled  there  by  replacement 
of  hot  primaries  with  cold  secondaries.  These 
processes  affect  the  discharge  stability  and  the 
thruster  efficiency. 

To  include  side  wall  effects  while  retaining 
the  simplicity  of  a  1-D  computation,  the  model 
described  here  starts  with  a  three-dimensional 
description  and  derives  a  1-D  approximation  in  a 
way  that  retains  the  side  wall  losses.  The  description 
is  time-varying,  includes  a  nonuniform  electron 
temperature  and  second  ionization,  uses  collision 
and  ionization  rates  that  depend  upon  Te  and  allows 
the  wall  surface  conductivity  to  limit  the  strength  of 
the  electric  field.  The  model  also  includes  ion 
pressure.  The  flow  of  plasma  into  the  walls  is  driven 
by  both  the  electron  and  ion  pressures,  so  both  are 
needed  to  calculate  the  wall  recombination. 
Although  one-dimensional,  the  model  allows  for  a 
diverging  channel  and  a  B  field  that  is  not  exactly 
radial.  It  calculates  the  ion  flux  into  the  walls  as 
well  as  the  incident  ion  angles  and  energies  as  a 
function  of  z,  which  can  be  used  to  predict  the 
erosion  of  the  walls  by  sputtering.  This  model 
assumes  quasineutrality,  ne  =  n  or,  including 
second  ionization, 

ne=n++2n++  [1] 

and  it  also  assumes  that  the  electron  gyroradius  is 
negligible  compared  to  other  distances  in  the 
problem.  These  two  assumptions  are  checked  during 
the  computation. 
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Since  the  curvature  in  a  thruster  channel  is  small 
relative  to  its  width,  local  Cartesian  coordinates  are 
used,  as  sketched  below. 


The  axial  (acceleration)  direction  is  “z”,  the  radial 
(across-channel)  coordinate  is  “x”  and  the  azimuth 
(with  which  nothing  varies)  is  The  x- 

dependence  is  approximated  analytically,  leaving 
equations  with  two  independent  variables,  z  and  t, 
that  are  solved  numerically. 

All  1-D  and  2-D  models  assume  azimuthal 
uniformity.  In  the  early  experiments  on  Hall 
thrusters,  this  was  far  from  true.  Probe 
measurements  showed  a  rotating  spike  of  density.1151 
Later  analysis  concluded  that  to  avoid  or  at  least 
limit  this  instability  one  needs  a  negative  axial 
gradient  of  (n/B).[16]  Models  like  this  one  assume 
that  this  is  done  correctly.  There  may  still  be 
azimuthal  waves  and  these  can  enhance  the  axial 
electron  transport.1'7,181  This  increased  conductivity 
is  included  in  the  model  via  the  anomalous  Hall 
parameter. 


2.0  Ion  and  Neutral  Dynamics 


The  ion  dynamics  are  described  by  the  Vlasov 
equations  for  singly  and  doubly  charged  xenon: 

dr  dr  eE,  dr  ,,  ■  ,  R  «.  „ 

4-  +  vr-  +  -7—  =  S{v-va)r^P ne  -  fne  Y 
at  M  ay 

[2] 


and 


dr  dr  2eEidT 
dt  drt  M  dvi 


=  fne7 


[3] 


Here/=/(x,y,t)  are  the  ion  distributions;  the  subscript 
i=x,  y  and  z;  a  repeated  index  means  a  sum;  ne  = 
n++2n++  is  the  electron  density,  nn  and  v0  are  the 
neutral  density  and  velocity  and  e  is  the  magnitude 
of  the  electronic  charge*.  The  right  sides  of  the 
equations  describe  ion  creation  and  loss  by 
collisional  ionization  of  neutral  xenon  and  second 
ionization  of  xenon  ions,  with  the  temperature- 
dependent  rates,  p  =  an<i  7~  7(Te>- 

Making  the  usual  definitions  of  density,  n,  velocity, 
it,,  pressure,  Pv  and  heat  flux,  Qm  in  terms  of  velocity 
moments  of  the  distribution  function,  one  obtains 
fluid-type  equations  by  taking  moments  of  [2]  and 
[3].  The  only  nonstandard  elements  are  the  source 
and  sink  terms. 

Simple  velocity  integrals  (“zeroth”  moments)  of  [2] 
and  [3]  give  continuity  equations, 

^-+  ^-(n+u*)  =  p  n0  ne  -  y  n+  ne  [4] 
dt  dri 

and  ^1+±(W+X++)  =  7„+„e  [5] 

dt  drt 

This  is  tensor  notation,  so  a  repeated  index  means  a 
sum  over  three  terms. 

Taking  the  first  moments,  and  using  the  continuity 
equations  to  simplify  them  gives  momentum 
equations, 


du *  + du, + 

n  M—T—  +  n  Mu,  1 


dt 


dr, 


dp; 


[6] 


eEjtl  '  dr, 


■  P M nen0(u  *  -v0j) 


and 


This  notation  could  be  confusing.  The  subscript  “z” 
doesn’t  mean  “ion”,  it  means  x,  y  or  z.  Ion  density, 
velocity,  etc.  are  denoted  by  superscripts:  n  ,  n++,  u+,  etc. 
The  ion  velocity  components  W/+  are  functions  of  position, 
but  the  components  of  in  equations  involving  /  are  not 
— -  they  are  the  other  three  coordinates  of  the  six¬ 
dimensional  phase  space.  Subscripts  are  used  for 
electrons  and  neutrals;  ne ,  n0 ,  ve,  v0.  Superscripts,  n  n  , 
etc.  would  be  more  consistent  with  the  ion  terminology, 
but  subscripts  are  more  conventional. 
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vTM- 


du ; 


d  t 


■  +  n*+Mu, 


dp,; 


du , 


dr, 


These  each  represent  three  equations,  j=  x,  y  and  z. 


[7] 


=  2  eE,n*+ - ^ - yMnen+(u* +  -«/) 

ari 

Taking  second  moments  of  [2]  and  [3]  and  using  [4]-[7]  to  simplify  gives 


and 


,d  +3  \  r>  +  , 

(37+m‘  +"^ 


+  K 

'Jk  dr 


iki  -  +  P*  ZZL.  +  pu+  7ZL.  +  P„ 


.  3  U  , 
dr. 


du : 


J‘  dr 

=  u0neM/3(Uj+  -v0j)(ut+  -vok) 


.3  ++  3  \  p  ++  .  ^Qjki 

(57+“'  fr)p“  “ST 


•+V 


3r.  w  dr  +  jl 


a«;+  _ 

dr 


These  each  represent  nine  equations. 


This  set  has  too  many  unknowns  to  solve.  The  usual 
simplification  is  to  disregard  the  Q^s  and  replace 
the  Py’s  by  scalar  pressures: 


=  7nep/  +  n*neM yiu*  -u*)(uk**  -u*) 

n0neM  J3(uj+  -  v0j)(uk+  -v0k)  =  0 
and  all  six  off-diagonal  equations  [9]  to 

rCneM  y(uj++  -Uj+)(uk++  -uk+)  =  0 


[8] 


[9] 


P+iJ=SvP+  P~*=S9P~  [10] 

Since  collisions  are  insufficient  to  isotropize  the 
pressures,  this  assumes  that  instabilities  do  so.  That 
may  not  be  true,  but  remember  the  objective,  to 
include  wall  losses.  The  ion  pressure  comes  from 
ions  being  made  at  different  places  and  therefore 
falling  through  different  voltage  drops,  creating  a 
velocity  spread  —  in  the  z  direction.  But  the  ion 
pressure  that  drives  the  flow  into  the  walls  is  in  the  x 
direction.  This  pressure  wouldn’t  even  exist  without 
collisions  to  isotropize  the  ions.  So  assuming  a 
scalar  pressure  is  the  worst  case,  where  the  ions  flow 
into  the  walls  as  fast  as  possible. 

Assuming  [10]  and  disregarding  Q  reduces  all  six 
off-diagonal  (zV j)  equations  [8]  to 


Since  uy+=  voy=  w/+=  0,  these  equations  are  trivially 
true  for  the  four  cases  where  either  j  or  k  is  y,  the 
azimuthal  coordinate.  The  remaining  two  cases, 
j=x;  k=z  and  j=z;  k=x  each  give 

n0neM  fKux+  -  vQx)(u2+  -  v0z)  =  0 

for  singly  ionized  and 

n+neM  y(ux++  -ux+)(u2++  -  uz+ )  =  0 

for  doubly  ionized  xenon.  Since  the  z  velocities  are 
not  zero,  these  equations  are  true  only  at  the  channel 
midplane,  where  the  x  velocities  are  zero. 
Everywhere  else,  the  best  one  can  say  is  that  the  x 
velocities  are  much  smaller  than  the  z  velocities,  and 
therefore  are  negligible. 


The  three  remaining  components  of  equation  [8],  namely  i=j-  x  and  i=j=  y  and  i=j-  z,  are 
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+  ^i_+  2p+  —f—  =  n0n,Mp  («/  -v0x)2 


and 


(bu:h;rn-)p''p  a,  a* 

(—  +  M,.+  ^-+  yne)p+  +  p+^  +  2p+  =  n,neMp  (u:+  -v0.)2 

3t  3r.  or, 


3z 


Adding  these  three  equations  gives 


3(—  +  u  +  —  +  yne)p*  +5p+^-=  n0neM/3  [(ux+  -v0r)2  +(w.+  -v0.)2] 
dt  or,  or, 

Using  continuity  [4]  to  eliminate  the  divergence  of  u+  gives 

^  8  5d+  „  *  +3  n+  8  n+ . 

i7~ir)= 

=  n0neMfl[(ux*  -v0x)2  +(u*  -v02)2] 


which  simplifies  to 


3  +  3  |  1 

—  +ui  ln 
,3t  3rJ 


(«+) 


5/3 


=  -^4{mh+[(u/ -v0J2  +(u..+  -v0:)2]  -  5p+}-rne 
3  p  n 

The  same  calculations  for  equation  [9]  give,  for  the  doubly-ionized  xenon: 


[11] 


3  + 

—  +  u, 

3 1  ' 


3  r, 


In 


u 


(O 


5/3 


n\r 

3p++n*+ 


Mri' 


+ 


3p+n+ 


[12] 


Equations  [11]  and  [12]  are  only  complicated 
because  of  the  source  terms.  Assuming  that  /?  and  y 
are  both  zero  collapses  them  to  p/(n5,i)—  constant, 
the  adiabatic  equation  of  state. 


even  include  ion  pressure.  It  is  included  here 
because  it  affects  the  flow  of  plasma  into  the  walls. 
However,  because  the  ion  pressure  is  a  small  effect, 
it  doesn’t  have  to  be  calculated  as  accurately  as  the 
other  variables. 


Altogether,  there  are  five  equations  for  each  species, 
a  continuity  equation  [4]  or  [5]  the  three  components 
of  the  momentum  ( F=Ma )  equation  [6]  or  [7]  and  a 
pressure  equation  [11]  or  [12].  The  latter  required 
some  approximations,  but  remember  that  the 
pressure  is  a  small  effect.  Most  1-D  models  don  t 


3.0  Quasi-One-D  Trial  Solution 

To  make  a  quasi-one-D  approximation,  we  assume 
that  nothing  varies  in  the  azimuthal  (y)  direction  and 
that  the  transverse  (x)  expansion  has  the  form, 
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Ux\x,z )  =  U*{x)-U*{z) 


[13] 


Simplified  versions  of  the  force  equations,  [6]  and 
[7]  are 


and  likewise  for  n+ ,  ux++,  p~,  etc.  To  obtain  a  trial 
form  for  ux+  (x)  etc.,  we  assume  that  the  plasma 

expands  uniformly  and  impinges  on  the  wall  at  the 
plasma  sound  speed,  vs.  This  is  consistent  with  the 
Bohm  sheath  criterion  that  plasma  enters  sheaths  at 
vs  or  faster.  Otherwise,  a  precursor  sheath  —  a 
rarefaction  fan  —  forms  to  accelerate  plasma  into 
the  sheath. 

The  problem  is  complicated  by  the  two  ion  species. 
The  usual  sound  speed  is  v5  s  -Jk(Te  +Ti)l  M  .  A 
generalization  is 


The  singly-charged  xenon  ions  reach  the  wall  at  vs+ 
while  the  doubly-charged  ions  reach  the  wall  at  v/+. 
This  is  an  approximation.  In  a  plasma  with  several 
species  of  ion,  there  are  several  sound  speeds,  but 
each  species’s  speed  involves  the  other  species. 
However,  these  forms  provide  an  overall  pressure 
balance: 


"=  ~k(Te+r)?f-  [16] 

dx  dx 

and 

:+r')^r-  t17l 

dx  dx 

Here  n  and  ux  are  functions  of  x  and  the  T  s  are 
constants.  Substituting  [13]— [15]  into  [16]  and  [17] 
gives, 

_  2x^ 

n+  ~  n0+e  w' 

and  [^] 

2xl 
++  w1 

n  =  n0  e  w 

So  in  the  x-momentum  equation,  a  linear 
acceleration  implies  a  Gaussian  density  profile. 
Since  the  walls  are  at  x  =  ±W/2,  the  density  at  the 
walls  is  e-0.5~  0.6  times  that  in  the  center  of  the 
channel.  Incorporating  the  no’s  into  the  z-dependent 
parts,  [18]  gives  the  trial  forms, 


/ j+M(02  +n+'M{V;+f  =  nekTe  +  p+  +  p++ 

and,  in  the  limit  of  negligible  ion  pressures,  they 
approach  v/+  =  -Jl  v/  as  they  should,  since  then 
the  ions  are  all  accelerated  through  the  same  E  field. 

Now  assume  that  both  species  expand  uniformly  in 
the  transverse  direction: 


lx1 

n+(x,z,t)=  n  +  (z,t)e  w’  [19] 

.rfl 

n++(x,z,t)=  n++(z,t)e  w'  [20] 

So  n(z)  is  the  density  at  x=  0,  at  the  midplane. 
Quasineutrality  gives  the  electron  density: 


where  W  is  the  channel  width.  In  these  coordinates, 
the  origin  is  at  the  center  of  the  channel  and  the 
walls  are  at  x  =  ±W/2.  In  general,  W  will  be  an 
increasing  function  of  z  because  the  channel 
diverges. 


ne(x,z,t)=  [«  +  (z,/)  +  2«++(z,o]  e  w  [21] 

Assuming  that  the  x-components  of  velocity  depend 
only  upon  x,  and  the  z-components  depend  only 
upon  z,  gives. 
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2v 

ix  (x,z,0-  ux  (x,0=  w  ■ X 

[22] 

N 

+ 

II 

N 

+ 

2v  ++ 

'■x  (X’Z’0-  U.r  (*>0=  w  X 

[23] 

L++(X,Z,0=  U:++(z,t) 

Finally,  for  the  pressures,  the  simplest  assumption  is 


p+(x,z,t)  =  n+(x,  z,t)kT*  (z,t) 
p++(x,z,t)  =  n++(x,z,t)kT++(z,t) 

So  the  pressure  is  obtained  from  density  and 
temperature  (i.e.  [24]  defines  the  ion  temperatures). 
The  density  is  described  by  [19]  or  [20]  and  the  ion 
temperature  is  assumed  to  vary  only  with  z,  not  x. 


4.0  Quasi-One-D  Derivation 

The  goal  now  is  to  substitute  these  forms  into  the  3-D  equations  to  get  a  one-D  description.  Logically,  Eqs. 

[19] _ [24]  are  guesses.  If  they  didn’t  lead  to  approximately  one-D  equations,  one  would  have  to  try  different 

forms.  But  as  will  be  seen,  very  little  further  approximation  is  necessary. 


Using  [19]  and  [22]  in  the  continuity  equation  [4]  gives 


t  a 

dt  dz 


J 


2v  + 
_ s _ 

W 


1- 


4xS 

W2 


[25] 


which  is  almost  x-independent.  Only  the  last  two  terms  still  vary  with  x,  and  that  variation  can  be  eliminated  by 
averaging  over  the  channel  or,  even  simpler,  by  considering  only  the  midplane,  x=  0.  Notice  that  if  one  just  wrote 
down  one-D  equations,  one  would  assume  a  continuity  equation  similar  to  this  -  except  for  the  last  term,  which  is 
new.  This  negative  term  is  an  additional  loss.  It  accounts  for  the  loss  of  plasma  to  the  side  walls  of  the  channel. 

Inserting  the  assumed  forms  into  equation  [5],  the  continuity  equation  for  the  doubly-iomzed  species  gives  a 
similar  result, 
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The  momentum  equations,  [6]  and  [7]  have  three  components,  but  the  y  components  are  trivial.  Assuming  that 
nothing  varies  withy,  these  just  give  0 =  0.  The  x  component  of  Eq.  [6]  gives 
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Neglecting  the  first  term,  which  depends  on  the  time  variation  of  the  temperature,  one  can  regard  the  rest  as 
specifying  the  ambipolar  E-field: 
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Since  we  are  only  interested  in  the  z-dependence,  we  don’t  really  need  this  equation,  except  to  check  for  accuracy 
of  the  model.  (If  Ex>  Ez,  then  the  problem  really  isn’t  quasi-one-dimensional.)  This  is  fortunate,  because  the  x- 
component  of  [7],  the  momentum  equation  for  the  doubly-ionized  species,  gives  a  slightly  different  equation  for 
the  ambipolar  E  field: 


Making  the  quasi-one-D  assumptions  in  the  z-component  of  [6]  gives, 

3?.* .  _  eE,  k  ar  tr  ar  p{r+i \ 

dt  2  dz  M  M  dz  Mn+  dz  n+ 

Except  for  a  possible  effect  of  n0,  this  has  no  x-variation  at  all. 

The  z-component  of  [7]  the  momentum  equation  for  the  doubly-ionized  species  is, 


[30] 


This  too  is  almost  completely  independent  of  x.  The  only  remaining  x  variation  is  in  the  last  term  (the  momentum 
change  due  to  production  of  Xe++  by  second  ionization  of  Xe+),  and  this  x-dependence  is  easily  eliminated  by 
looking  at  the  center  of  the  channel.  Again,  the  result  is  not  very  different  from  what  one  would  get  by  merely 
writing  down  a  one-D  description,  except  that  it  now  includes  a  drag  due  to  ion  creation. 

So  far,  the  quasi-static  approximation  has  been  surprisingly  successful.  The  assumption  of  a  uniform  expansion 
in  x,  plus  the  Gaussian  density  profiles  that  match  such  an  expansion  has  essentially  eliminated  the  x-dependence, 
leaving  equations  that  depend  only  on  z  and  t.  All  that  remains  is  to  approximate  the  pressure  equations. 

Eq.  [1 1],  the  pressure  equation  for  the  singly-ionized  xenon  gives, 
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the  x-dependence  has  certainly  not  disappeared  here,  but  there  is  no  reason  why  it  should.  The  assumed  x- 
variations  of  velocity  and  density  were  chosen  to  simplify  the  continuity  and  momentum  equations,  not  this 
equation.  But  notice  that  the  left  hand  side  is  independent  of  x  and  the  right  side  can  be  greatly  simplified  by 
looking  only  at  x=0. 
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Inserting  the  quasi-one-dimensional  forms  into  [12],  the  pressure  equation  for  doubly-ionized  xenon  gives, 
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This  equation  is  similar  to  [31]  but  simpler,  because  there  is  no  collisional  loss  of  these  ions.  It  also  simplifies  at 
x=0.  This  completes  the  model  of  the  ion  dynamics. 


Restating  the  whole  set  of  ion  equations  for  flow  along  the  midplane  of  the  channel,  x-0 
depend  only  on  z  and  t,  as  intended.  From  [25],  the  continuity  equation  for  Xe+,  we have 
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And  [26]  at  x=0  gives  the  continuity  equation  for  Xe*+, 
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The  momentum  equations  were  essentially  independent  of  x  as  obtained  already.  From  [29]  we  have,  for  the 
singly-ionized  xenon, 
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[35] 


The  only  change  from  [29]  is  that  the  neutral  density  and  velocity  have  been  decorated  with  double  tildes  to 
signify  that  any  x- variation  in  them  has  been  eliminated  by  taking  the  values  at  the  midplane  of  the  channel. 

Equation  [30]  taken  at  x=0  gives  the  one-D  momentum  equation  for  the  doubly-ionized  xenon, 
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From  [31]  we  have  the  equation  for  the  pressure  —  i.e.  for  the  temperature  —  of  the  singly-ionized  xenon  at  the 
midplane, 
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Finally,  [32]  at  x=0  gives  a  similar  equation  for  the  temperature  of  the  Xe++, 
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These  six  equations  give  the  ion  dynamics.  Since  there  are  six  variables,  namely  the  densities,  z-velocities  and 
temperatures  of  the  two  species,  six  equations  are  enough.  These  equations  do  not  explicitly  include  energy 
equations  for  the  ions,  but  these  are  implied  and  can  be  derived  from  the  foregoing  results.  But  because  the 
energy  equations  are  more  complicated,  it  makes  sense  to  use  the  above  forms. 

Remember  that  these  final  results  are  all  at  the  midplane  of  the  channel,  in  the  quasi-one-dimensional 
approximation,  [19]  —  [24],  So  the  total  flux  per  unit  circumference  is  not  nuzW,  but 
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[39] 


The  density  at  the  midplane  is  n  ,  the  density  at  the  walls  was  shown  by  Eq.  [18]  to  be  0.6  n  and  Eq.  [39]  shows 
that  the  mean  density  averaged  across  the  channel  is  0.86«  . 


This  completes  the  model  of  the  ion  dynamics.  The  equations  that  describe  the  other  species  in  the  plasma,  the 
electrons  and  the  neutrals  will  be  considered  only  in  their  one-dimensional  forms  along  the  midplane.  Since  all 
the  quantities  in  the  following  equations  refer  to  the  midplane,  the  double  tildes  will  be  omitted  from  here  on. 


5.0  Neutrals 


The  equation  of  continuity  for  the  neutral  gas  is 


3«0  3 

3  t  3 z 


(«0V0z) 


=  -  0no  ( n +  +  2 «++)  + 


2(v,y +vro 


[40] 


W 


The  last  term  describes  neutrals  coming  off  the  wall 
where  plasma  has  recombined.  For  consistency,  the 
neutral  density  distribution  in  x  is  assumed  to  be 
Gaussian  too.  This  conserves  particles  both  along 
the  centerline  and  across  the  channel.  But  it  is  an 
approximation.  Because  neutrals  are  coming  off  the 
walls,  and  because  the  rate  of  loss  by  ionization  is 
greatest  in  the  center  where  the  plasma  density  is 


highest,  the  neutral  density  really  peaks  at  the  edges 
of  the  channel,  but  modeling  that  correctly  would 
require  a  full  two-D  analysis.  Putting  too  many 
neutrals  in  the  center,  which  this  does,  makes  the 
ionization  rate  too  high.  To  compensate  for  this,  the 
ionization  cross  section  could  be  reduced. 


Since  neutral-ion  collisions  are  rare,  the  neutrals 
diffuse  along  the  channel,  colliding  only  with  the 

1  3  n0  1 

walls.  For  a  small  density  gradient, - —  s  — , 
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and  the  neutral  gas  temperature  T0  is  constant. 


For  sharp  density  gradients,  the  velocity  defined  by 
[41]  becomes  much  greater  than  the  thermal 
velocity,  which  is  unrealistic.  To  avoid  this,  the 
neutral  velocity  will  be  estimated  as 
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Hall  parameter,  HP,  and  writing  the  electron  velocity 
as 
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The  electrons  also  satisfy  a  continuity  equation, 
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Since  the  electron  density  is  determined  by 
quasineutrality,  this  is  a  constraint  on  the  electron 
velocity,  and  therefore,  through  Eq.  [44],  on  the 
electric  field. 


7.0  Collisions 


6.0  Electron  Dynamics 

The  electron  density  is  determined  by 
quasineutrality,  Eq.  [1].  The  electron  axial  velocity 
is  the  collisional  drift  across  the  B-field  due  to  the 
electric  field  and  the  pressure  gradient: 

(I classical )  1  77  ^  d(tlekTe) 

v  = - t.  H - - - 

B  [  ene  dz 

where  Vj  and  v2  are  elastic  and  inelastic  collision 
frequencies  and  £2e  =  —  is  the  electron 

me 

gyro  frequency.  The  ratio  (V/+  Vi)/Qe  is  the 
reciprocal  of  the  Hall  parameter.  This  is  the  same 
equation  used  in  other  models/"'  except  that  here  it 
is  in  SI  units. 

The  Hall  parameter  may  be  either  classical  or 
anomalous.  Considering  Vj  and  v2  to  be  the  classical 
collision  frequencies,  one  knows  that  Hall  thrusters 
can  operate  in  regimes  where  this  would  give  an 
unrealistically  high  Hall  parameter.  In  that  case, 
plasma  fluctuations  (e.g.  azimuthal  waves)  that 
scatter  electrons  across  the  magnetic  field  will  be  the 
primary  cause  of  electron  current.  This  can  be 
included  in  the  theory  by  assuming  an  anomalous 


Vi+Vi 


The  ionization  rate  /?,  the  second  ionization  rate  y 
and  the  electron  elastic  (V;)  and  inelastic  (v2) 
collision  frequencies  are  all  functions  of  the  electron 
temperature  Te.  Fortunately,  some  data  are  available 
for  the  energy  range  of  interest/191  The  simplest 
procedure  is  to  approximate  these  data  with 
empirical  formulas. 

The  elastic  collision  frequency,  vh  is  the  sum  of  the 
frequency  of  collisions  with  neutral  atoms,  singly- 
charged  ions,  and  with  doubly-charged  ions, 

v,  =  Rlon0  +  n++  [46] 

R,o  is  the  rate  for  electron  collisions  with  neutrals 
and  Rf+  and  Ri++  are  the  Coulomb  collision  rates. 

The  electron  inelastic  collision  frequency  also 
includes  collisions  with  neutrals  and  with  singly  and 
doubly  charged  ions': 

_  -^20w0  +  ^  ^2  W  J-47J 

2  kTe 

So  one  needs  either  tables  or  formulas  for  eight 
functions,  namely  the  two  ionization  rates,  P(Te)  and 
y 'TJ  (in  m3/s),  the  three  elastic  collision  rates, 

RnifTJ  R,+(Te)  and  R^fTe)  (also  in  ms/s),  and  the 
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three  inelastic  rates,  Rio(Te)  Ri  (T^z. nd  i?2++(T&)  (in 
W-m2). 


Solving  these  equations  gives  the  sheath  voltage: 


8.0  Side  Wall  Boundary  Conditions 

If  the  walls  are  not  conductive,  the  net  current  into 
the  walls  must  vanish,  so  the  electron  and  ion  wall 
fluxes,  Te,  r  and  T'f+  have  to  balance: 

re  =  r+2r++  +  sre 

where  S  =  S(Te)  is  the  wall’s  electron  secondary 
emission  coefficient.  Since  the  electron  thermal 
velocity  is  higher  than  that  of  the  ions,  the  electron 
flux  would  be  greater,  unless  balanced  by  secondary 
emission  or  reduced  by  a  sheath  voltage  barrier.  The 
ion  fluxes  are 

r=o.6«v  r^o.bn+v* 

The  factors  of  0.6  enter  because  the  densities  at  the 
channel  edge  are  only  0.6  times  those  in  the  center 
where  n+  and  n  are  defined.  According  to  the 
analysis  of  the  MIT  group'61  the  primary  electron 
flux  into  the  wall  is 

fgv 

T  =  0.6(«++2«++)^e  kT‘  [48] 
4 

where  <Pw  is  the  sheath  voltage  (negative  because  the 
electric  field  points  toward  the  wall),  ce  is  the  mean 
electron  thermal  speed'201: 

[49] 


and  me  is  the  electron  mass. 

The  effective  secondary  emission  coefficient,  Sejj  is 
derived  in  Ref.  6  to  be 


8=  T(2  +  B)A0 


0.198 


[50] 


where  Ao=0.141  and  B0=0.576  for  boron  nitride  and 
IJx)  is  the  Gamma  function. 
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If  the  ion  temperatures  are  negligible,  v/-v/+  and 
Eq.  [61]  reduces  to 


«,  =  «LJ 

'  e  U(l-£) 


Given  the  two  densities  and  the  three  temperatures, 
Eq.  [51]  predicts  the  sheath  voltage.  Because  ce» 
v/,  vs++,  the  argument  of  the  logarithm  is  much  less 
than  one  and  the  sheath  voltage  is  negative  —  unless 
the  secondary  coefficient  Jis  nearly  one.  A  negative 
voltage  means  the  E  field  points  toward  the  wall. 

Eq.  [50]  says  that  8  goes  to  one  as  the  temperature 
approaches  1.93xl05  °K,  which  is  16.64  eV.  At  a 
slightly  lower  temperature  the  argument  of  the 
logarithm  in  Eq.  [51]  becomes  one,  driving  the 
sheath  voltage,  <ZV  to  zero.  One  does  not  expect  Te 
to  get  this  high  because  secondary  emission  cools 
the  electrons  by  replacing  hot  primaries  with  cold 
secondaries.  Indeed,  this  entire  description  is  valid 
only  if  the  sheath  voltage  &w  is  negative. 

Ions  strike  the  wall  with  the  kinetic  energy  they  had 
upon  reaching  the  sheath,  plus  the  energy  gained  by 
falling  through  the  sheath.  For  singly-ionized 
xenon,  this  energy  is, 

[53J 

and  doubly-ionized  xenons  strike  the  wall  with 
energy, 

V~=  1a /(v/+)2  +2e*w  +$At(u~f  [54] 

These  energies  and  the  fluxes  I*  and  V  can  be 
used  to  estimate  the  erosion  of  the  wall  by 
sputtering.  But  this  picture  assumes  that  the 
acceleration  of  ions  by  the  discharge  voltage  is 
precisely  tangential  to  the  wall,  which  is  true  only  if 
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the  5  field  is  exactly  radial.  In  reality,  the  5-field 
lines  are  curved,  so  even  if  the  discharge  E  is  normal 
to  5,  it  can  drive  ions  toward  or  away  from  the 
walls.  But  as  discussed  below,  the  model  can  be 
corrected  for  this  effect. 


K=  1.7-1  O'13  r//2  W/cm-deg.  For  the  conditions 

in  a  Hall  thruster,  this  keeps  Te  quite  uniform  across 
the  channel.  The  electron  temperature  is  therefore 
determined  by  the  energy  balance,  integrated  across 
the  channel 


9.0  Electron  Temperature 

The  only  piece  still  missing  is  the  calculation  of  Te. 
This  must  be  obtained  from  the  electron  energy 
balance,  the  requirement  that  the  change  in  the 
electron  heat  content  equals  the  heating  minus  the 
cooling.  In  fact,  the  heat  content  is  negligible,  so  the 
electron  temperature  must  be  whatever  makes  the 
heating  and  cooling  almost  equal,  but  for 
computation,  it  is  simpler  to  use  the  energy  equation. 
The  heating  equals  the  electric  field  times  the 
electron  current.  Cooling  is  by  inelastic  collisions 
and  by  heat  loss  to  the  walls,  where  hot  primaries  are 
replaced  by  cold  secondaries. 
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where  qin  is  the  electrical  power  input  to  the 
electrons,  qcoa  is  the  power  lost  by  inelastic 
collisions  with  atoms  and  ions,  both  integrated 
across  the  channel,  and  q^n  is  the  power  lost  to  the 
walls.  All  three  of  these  q’s  vary  with  the  electron 
temperature. 

Integrating  across  the  channel,  as  in  Eq.  [39]  gives 
the  power  input: 

q.m  =  -  0.86  W(n+  +  2 n*+)evezE  [56] 


For  wall  cooling  to  affect  Te  near  channel  center,  the 
heat  must  conduct  through  the  plasma.  For  thruster 
conditions,  Spitzer’s  result[21)  for  the  thermal 
conductivity  of  the  plasma  gives 

The  collisional  power  loss  is  quadratic  in  the  density: 


Here  va  is  the  electron  cross-field  drift  velocity 
given  by  Eq.  [44]  which  does  depend  on  Te. 
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[57] 


where  V2  is  the  inelastic  collision  rate. 

The  electron  energy  loss  to  the  walls  is  shown  in  Ref.  6  to  be 

2T,[lKT,-ST^)-i\-S)e<S>A  l58) 

where  Tsec  is  the  temperature  of  the  emitted  secondary  electrons  (approximately  1  eV)  and  a  factor  of  two  has 
been  added  to  the  result  in  Ref.  6  to  account  for  the  loss  to  both  the  inner  and  outer  walls.  Note  that  qin,  qcou,  and 
qwaii  are  all  in  W/m2. 

Using  these  in  Eq.  [55]  gives  an  equation  for  electron  temperature: 
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Besides  the  explicit  Te  dependence  here,  the  drift  velocity,  ve  (Eq.  [44]),  the  thermal  velocity  ce  (Eq.  [50])  the  wall 
sheath  potential,  &w  (Eq.  [52])  and  the  secondary  emission  coefficient,  8 (Eq.  [50])  are  all  functions  of  Te. 

For  a  numerical  solution,  this  can  be  simplified,  because  the  electron  heat  capacity  is  small  compared  to  the 
heating  and  cooling  rates.  Physically,  Te  assumes  the  value  that  makes  the  net  heating  (the  rhs  of  this  equation) 
vanish.  Eq.  [59]  is  simply  a  way  to  find  that  temperature.  To  avoid  introducing  shorter  times  and  higher 
velocities  into  the  calculation,  the  heat  convection  was  omitted  and  the  heat  capacity  artificially  increased  in  the 
computations. 


10.0  Non-Axial  Flow,  Non-Constant  Width 

With  the  foregoing  equation  for  electron 
temperature,  the  model  is  complete.  However,  as 
now  formulated  it  contains  an  inaccuracy  that  would 
limit  its  utility,  namely  the  assumption  that  the  walls 
are  straight  and  B  is  normal  to  the  walls.  Where  B- 
field  lines  bow  outward,  an  E-field  normal  to  B  can 
point  slightly  toward  the  wall  and  ions  moving 
normal  to  B  can  strike  the  wall,  causing  additional 
loss  of  plasma  and  additional  sputtering  of  walls. 
This  is  worth  including. 

A  diverging  B  also  changes  the  cross-sectional  area 
of  a  flux  tube  from  wall  to  channel  center.  The  wall 
losses  (the  0. 6  n+vs+  used  above)  must  be  corrected 
for  this  area  ratio.  Because  B  is  stronger  at  the  wall, 
the  cross  sections  of  flux  tubes  are  smaller  there  and 
therefore  the  loss  of  ions  to  the  wall  is 
proportionately  less. 

As  shown  in  the  drawing  below,  one  can  define  an 
angle  y/thsX  a  normal  to  B  makes  with  the  wall  and 
calculate  the  additional  wall  loss  that  this  causes. 
(For  simplicity,  it  is  assumed  that  a  field  line  makes 
the  same  angle  with  the  inner  and  outer  walls.)  One 
can  also  allow  the  channel  width  to  vary  —  that  is, 
W  becomes  W(z)  —  and  correct  the  plasma  density 
for  the  additional  expansion.  The  reduction  in  the 
wall  losses  due  to  the  area  being  smaller  at  the  wall 
can  be  included  through  an  area  ratio,  A(z)^(area  at 
walls)/(area  at  midplane),  which  in  general  will  be 
less  than  one.  Finally,  one  can  recognize  that  the  B- 
field  strength  varies  with  z.  That  is,  to  solve  the 


model,  first  input  four  functions,  y/(z),  W(z),  A(z) 
and  B(z)  that  describe  the  geometry  of  the  channel 
and  the  field.  The  convention  is  that  z  denotes  the 
distance  along  the  midplane  of  the  channel. 
Distances  z  along  walls  are  not  physical  distances, 
but  mappings  of  the  corresponding  midplane 
distances,  with  the  mapping  generated  by  following 
5-field  lines. 

If  the  B  field  curves  as  sketched,  ion  trajectories  will 
not  necessarily  be  exactly  normal  to  B,  even  if  the  E 
field  is.  However,  the  E  field  generally  increases 
with  z,  so  the  ions’  direction  is  primarily  affected  by 
the  E  field  near  their  location,  which  is  normal  to  B. 
So  we  assume  that  the  ion  “z-velocity”  is  exactly 
normal  to  B,  i.e.  directed  into  the  wall  at  an  angle  y/ 
from  tangential. 


Figure  2  Definition  of  the  angle  v|/. 


If  this  angle  were  negative,  as  it  is  in  the  lowest  field 
line  drawn  above,  the  ion  velocity  would  be  angled 
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away  from  the  wall.  This  can  obviously  happen,  and 
would  be  beneficial,  because  it  would  move  ions 
away  from  the  wall,  but  is  also  outside  the  scope  of 
this  model,  which  assumes  the  plasma  has  a 
Gaussian  profile  (Eqs.  [19]— [21]).  In  this  formalism, 
a  negative  ys  would  make  the  wall  a  source  of 


plasma,  which  is  not  physical.  Therefore  negative 
values  of  \tfz)  will  not  be  input.  If  y/  is  really 
negative,  as  on  the  lowest  field  line  in  the  above 
drawing,  it  will  be  set  to  zero.  Note  also  that  y/  is 
defined  with  respect  to  the  wall,  not  to  the  vertical. 


If  the  velocity  normal  to  B  is  v.  and  that  parallel  to  B  is  vs,  then  the  total  velocity  is 
vT  -  -Jv}  +  v/  .  As  illustrated  to  the  left,  this  total  velocity  makes  an  angle 
£  =  tan'l(vi  /v.)  with  vr.  The  angle  between  v-  and  the  wall  is  ( yA^Ej  and  the  velocity 

normal  to  the  wall  is  v±  =  vT  sin(^  +  £) .  The  flux  of  ions  into  the  wall  is  7V=  v±n.  For 
sputtering,  this  flux  (for  each  species)  strikes  the  wall  at  velocity  vy,  at  an  angle  ( yfr  q). 
But  this  is  not  the  right  wall  loss  to  use  in  the  continuity  equations,  because  7" v  is  the  flux 
per  unit  area  of  wall.  For  continuity,  one  needs  the  loss  per  unit  area  normal  to  B.  That 
flux  is  F=  rw/cos( y/),  which  replaces  vsn  in  the  continuity  equations. 


The  other  needed  modification  is  the  effect  on 
density  of  a  diverging  channel,  which,  to  conserve 
particles,  must  obey 


needs  (l/n)(dn/dt)  =  (flux)/(total  particles)  = 
1.40(vl/W).  So  the  coefficient  of  the  last  term  in 
the  continuity  equation  should  be  1.4  and  not  2. 
Then  the  continuity  equation  becomes. 


l  dn_  1  dW_ 
n  dz  W  dz 

Making  these  changes  in  the  ion  continuity  equation 
[33]  gives 


+  sin(£+  +  yf)  I  +2  +2 

where  v,  = - -Ju,  +v, 

cos  (y/) 

f..+  \ 


and  =  tan 


\Uz  j 


dn*  1  9  I ,,,  +  +i 
dt  WdzX  ■  ' 


(j3n0-rn'){n*+2n") 


L4  AvL+n+ 
W  "" 
[60] 


For  calculation  of  sputtering,  the  flux  per  unit  area  of 
wall  is  = 0.6  n+ v i  cos  ( y/) 

For  the  doubly-ionized  xenon,  [34],  the  same 
calculation  gives 


dn 


+  ——[wn”u”)=  yn*{n*  +  2n^)  - 
dt  WdzK  '  ' 


XAAv~n* 

W 

[61] 


This  continuity  equation  was  derived  by  making  the 
quasi  one-dimensional  approximation  and  then 
setting  x=0.  One  way  in  which  this  is  approximate 
is  that  it  does  not  exactly  conserve  particles.  To  see 
this,  note  that  the  total  flux  of  particles  into  the  two 
side  walls  of  the  channel  is  2(vL+)(0.6  n). 
Integrating  across  the  channel  gives  the  total  number 
of  particles,  0.86(Wn).  To  conserve  particles,  one 


where  v~  =  sin(<T_t ~)2  +(v/+)2  and 
cos(^) 

f  ++^ 

£++  s  tan-1 


++ 

\u--  j 


and  A=A(z)  is  the  area  ratio  defined  above. 
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those  equations  will  change  because  of  the  changes 
The  momentum  [35]&[36]  and  pressure  [37]&[38]  in  the  densities  caused  by  the  extra  terms  in  Eqs. 

equations  are  not  changed  in  form  by  these  effects,  [60]&[61]. 

although  the  velocities  and  pressures  predicted  by 


Making  the  same  additions  to  the  neutral  continuity  equation  [40]  gives 


N  *  ,  ♦  .  „  .  1  AA(yL*n++vL"n*) 

1 f+F&(’rv’*)='M("  )+ - w - 


[62] 


Finally,  the  electron  continuity  equation,  [45]  including  these  effects  becomes 


l^L  +  ^T('Wn'V^=  (^no+7n+)ne  ~^-(vL+n+  +2vi++n++) 
dt  W  oz  W 


[63] 


11.0  Summary  of  the  Model 

This  completes  the  formulation.  There  are  twelve 
variables,  namely  the  four  densities  and  the  four  z- 
velocities  for  the  four  species  —  the  neutral,  singly- 
and  doubly  charged  xenon  plus  the  electrons  —  the 
two  ion  temperatures,  the  electron  temperature,  and 
the  electric  field  normal  to  B.  (The  neutral 
temperature  is  a  constant,  not  a  variable.)  These 
twelve  variables  are  defined  in  the  two-dimensional 
(z,t)  space. 

There  are  also  twelve  equations,  quasineutrality,  Eq. 
[1],  the  continuity  [60]&[61],  momentum  [35]&[36] 
and  pressure  [37]&[38]  equations  for  both  ion 
species,  the  continuity  [62]  and  diffusion  [43] 
equations  for  the  neutrals,  Eq.  [44]  for  the  electron 
drift  velocity,  Eq.  [63]  for  electron  continuity  and 
Eq.  [59]  that  gives  the  electron  temperature 
implicitly.  These  equations  in  turn  use  the 
collisional  rates,  the  electron  secondary  emission 
coefficient  of  the  wall  material,  Eq.  [50]  and  the  wall 
sheath  voltage  drop,  Eq.  [52]. 

But  it  is  not  necessary  to  solve  this  many  equations 
numerically.  The  quasineutrality  condition,  Eq.  [1] 
gives  the  electron  density  in  terms  of  the  ion 

densities:  ne  =  n+  +  2n++ ,  which  can  be  used  to 
eliminate  ne  from  the  other  equations. 


—  [w(n+uz+  +ln++u2++  -(n+  +  2«++)vfZ)]=  0 
az 

which  says  that  the  axial  current  is  conserved.  That 
is, 

n*u*  +  2n**u+*  ~(n+  +2 n++)ve: 


where  J(t)  is  the  total  discharge  current  per  unit 
circumference  of  the  annulus. 

This  current,  J(t),  is  an  input  to  the  model,  but,  as 
discussed  below,  one  can  model  an  external  circuit 
that  is  essentially  a  voltage  sulrce  by  an  inversion 
procedure.  Once  J  is  obtained,  this  equation,  which 
replaces  the  electron  continuity  equation,  [63]  can  be 
solved  for  the  electron  drift  velocity: 


+  +./■»++  ++ 
n  u,  +2 n  u,  - 


J 

eW 


n+  +2n 


[64] 


and  this  can  be  used  to  eliminate  from  all  the 
other  equations. 


The  electron  and  ion  continuity  equations,  together 
with  quasineutrality,  imply 


15 


Equation  [44]  for  the  electron  drift  velocity  can  then 
be  solved  for  the  electric  field: 


E=  -  5vc 


1-1/2 


V  J 


(H„y 


1  d(jtJcT\)  |-g5j 
ene  dz 


Finally,  adding  the  continuity  equations  [60],  [61]  and  [62],  gives 

-(fl++«++  +  K0)  +  —  |~M«V  +  «+V+ +«<,%)]  =  0  [66] 

dr  wdzv  v  -  2  ' 

which  shows  that  the  flux  of  particles  through  the  thruster  is  conserved.  This  can  be  used  instead  of  Eq.  [62],  the 
continuity  equation  for  the  neutrals. 

With  these  simplifications,  the  model  consists  of  eight  variables,  n+,  uE  and  T ;  n+,  uz++  and  7rf+;  n0  and  Te.  The 
electron  density  and  drift  velocity,  the  neutral  velocity  and  the  electric  field  are  functions  of  these  primary 
variables.  The  eight  primary  variables  obey  eight  equations,  namely  [60],  [35]  &  [37]  for  Xe+;  [61],  [36]  and  [38] 
for  Xe++;  [66]  for  the  neutrals  and  [59]  for  Te.  One  also  needs  the  functional  forms  of  the  collision  rates  plus  8 
and  &w,  and,  finally,  a  model  of  the  external  circuit. 

Since  the  variables  are  functions  of  z  and  t,  the  evident  procedure  is  to  store  all  the  needed  quantities  as  functions 
of  z  (e.g.  at  1000  points)  and  advance  the  system  in  time.  Besides  the  eight  primary  variables  and  four  secondary 
variables  (ne,  va,  vfe  and  E),  there  are  also  the  input  functions  y/(z),  B(z),  A(z)  and  W(z)  that  describe  the  channel 
and  5-field  geometry. 

While  solving  these  equations,  one  wants  to  check  that  the  two  conditions  mentioned  earlier  are  satisfied,  namely 
(1)  that  the  electric  field  and  electron  density  are  consistent  with  quasineutrality  and  (2)  that  the  electron 
gyroradius  is  small  compared  to  distances  over  which  the  electron  density  changes.  This  is  conveniently  done  by 
defining  two  ratios, 


ene  dz 


5,  = 


dn. 


£2.  n _  dz 


Both  of  these  should  be  <1  or  the  model’s  assumptions  are  not  satisfied. 


[67] 


12.0  Reformulation  of  the  Equations 

The  next  step  is  to  put  the  equations  into  a  form  suitable  for  numerical  solution.  To  integrate  in  time,  one  stores 
eight  variables  -  the  density,  z-velocity  and  temperature  of  each  ion  species,  plus  the  neutral  density  and  the 
electron  temperature  -  as  functions  of  z.  From  these  eight,  the  other  four  variables  -  vfej  E,  ne,  and  va  -  can  be 
calculated,  as  well  as  all  the  z-derivatives.  Then  one  can  calculate  the  time  derivatives  of  the  eight  primary 
variables,  and  advance  the  system  in  time. 

Restating  the  equations  for  these  eight  primary  variables,  one  has,  from  Eq.  [60], 

—  =  Ft  [68] 

dt 


where 
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=  -^{Wn+uz+)  +  (J3no-r"+)(n+  +  2«++)  -1AA£-— 


And  Eq.  [35]  gives 


3  us* 
3 1 


1  2 
M 


where 


,  M/3(rt+  +2n++)n0(uz+  -v0z) 


-|-(00 

n„  oz 


and 


F„+»  |a/(iii+)2+*3T*+^ 


For  T,  Eq.  [37]  gives 


1  3  r  2  3m + 


r+  3t  3m +  3? 


=  F 


where 


F,  =  -M.+  ^-ln 
3z 


f  kT*  N 
(«+)2/: 


Mn0P 


3kT* 


1  + 


2m4 


A 


(«-+-v0.-)2- 


2  5  kT 


M 


-Yin*  +2n++) 


For  m++,  Eq.  [61]  gives 


a  n++  _ 


where 


F4h  — (jFm++m.")+?'«+(«+  +  2'j++) 

4  WdzK  '  ’  V  7  fF 


For  the  velocity,  from  Eq.  [36], 


3m. + 


3? 


_2_5_ 

M 


where 


F5  =  — (O)  +  *r+  9lD^---  -  20  +  2kTe 

3z  or  3z 

M;Kn++2n+»z++-0 


2d 

—  —{kTene) 
• m_  dz 


and 


F„++  =  2-M  (uz++)2  +kT*  +2kTe 


For  O  Eq.  [38]  gives 


i  3 r+ 2_^!1=  f 

r++  3/  3m  ^  3f  6 


[69] 

[70] 

[71] 

[72] 

[73] 

[74] 

[75] 

[76] 

[77] 

[78] 


17 


where 


F‘"-M=  3z 


3  ,  (  kT++  )  .  M(n+  +2n++)n+r 


++  w  In 


(«++)2 


+  - 


3  kT++n* 


(«.  -M.  )  “ 


+,2  5/tr+  .  3/tr 


M  M 


For  the  neutral  density  we  have,  from  Eq.  [66] 


K.  +  ^l  +  ^  =  _F 

3t  3?  3 t 


where 


f7=  -^-^;{w(n+u*  +n++u*+  +n0V0z)] 


Finally,  for  Te,  we  have,  from  Eq.  [59],  — 

n _  at 


where 


■nkT 

e  e 


=  -F, 


vK  3  f  3  ^ 

F„  =  — - —nkT 

8  3zl  2  e  e 


+  evezE  +  0.87v2kTe 


A  s* 

+  0.35  -f-  e  kT‘  [2k(Te  -  STsk )  - (1  -  $)e<bw ]  = 


v.  3  fl  A 
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-nkT. 


-  eBvi 
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+- 
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/  H ; 


+  0.87v2kTe 


+  ■ 


1.4  A 
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n*v*  +  2n*+vs 

(1  -S)n, 


[2k{Te-ST^)-{\-S)e<S>w] 


For  numerical  solution,  it  is  convenient  to  unscramble  the  mixed  equations  in  this  set,  writing 

=  rF  |  2t+fx 

3 1  3  3 n+ 


instead  of  Eq.  [72]  and 


MZ  -  r+  fa  i  2T++F< 
3 1 


6  ++ 
3  n 


instead  of  Eq.  [78], 
instead  of  Eq.  [80],  and 


£^  = _ L. _ <F  +2F  )  -  35- 

3 1  n+  +  2n*+K  '  *’  3 k 


instead  of  Eq.  [82]. 


[79] 

[80] 

[81] 

[82] 
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To  advance  the  system  in  time,  starting  with  the  eight  primary  variables  as  functions  of  z,  and  therefore  then- 
gradients  in  z,  one  must  calculate  the  neutral  velocity  from  Eq.  [43]  and  the  collisional  rates  from  the  electron 
temperature.  Then  one  can  calculate  the  eight  F’s  and  therefore  the  eight  time  derivatives  needed  to  advance  the 
system. 

In  practice,  these  equations  as  they  stand  give  difficulties,  because  they  predict  sharp  gradients  that  cause 
numerical  problems.  This  is  not  physical  because  particle  diffusion  smooths  the  variables.  So  it  makes  sense  to 
add  diffusion  to  the  model.  Taking  a  mean  free  path  of  the  order  of  the  channel  width,  while  limiting  the 
diffusion  velocities  to  realistic  values  gives  diffusion  coefficients, 


D+  = 


kT+ 

r  i 

dn+ 

3  ' 

V  M 

l n+ 

dz 

+  w) 

D*+  = 


lfcT++ 

(  1 

a«++ 

3  ] 

4- 

1  M 

l"++ 

dz 
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In  the  calculation,  these  were  used  for  the  densities,  while  the  coefficients  for  the  velocities  and  temperatures  were 
taken  as  one  half  and  one  fifth  of  these  rates,  respectively. 


Mathematically,  these  additions  are  an  estimate  of  the  higher  terms  in  the  velocity-moment  expansion  that  were 
neglected  in  the  derivation.  But  to  derive  them  that  way,  one  would  have  to  remember  that  the  system  does  have 
walls,  ion-ion  collisions  and  fluctuating  electric  fields  in  the  plasma.  Practically,  the  effect  of  the  diffusion,  aside 
from  making  the  calculation  stable,  is  small.  The  diffusion  terms  in  one  run  were  increased  by  50%.  The 
discharge  current  did  increase,  but  only  by  9%. 


13.0  External  Circuit 

The  thruster  driver  was  modeled  as  a  voltage  source  with  a  series  resistance  and  a  parallel  capacitance.  The 
discharge  current  (which  the  model  uses)  and  the  discharge  voltage  (which  the  model  calculates)  imply  the  current 
through  the  capacitor  and  therefore  the  derivative  of  the  voltage  on  the  capacitor  and  discharge.  Since  the 
discharge  current  is  input  to  the  model,  one  needs  the  change  in  current  that  will  give  the  voltage  change  required 
by  the  circuit.  In  calculating  that  needed  dl/dt,  one  must  remember  that  the  discharge  resistance  varies  with  time. 
The  needed  equation  is 


dV__  dV_  dl_  |  dV 
dt  dl  ,  dt  dt 


or, 


dl_ 

dt 


dV_dV 
dt  dt 
dV\ 


[88] 


To  evaluate  this,  the  program  calculates  three  E  fields  and  resulting  voltages:  the  V due  to  the  actual  current,  the  V 
due  to  a  slightly  different  current  (for  dV/dI\,)  and  the  V  for  the  current  on  the  preceding  time  step  (for  dV/dt\i).  In 
practice,  the  discharge  voltage  holds  constant  to  within  a  fraction  of  a  percent,  so  the  power  supply  is  effectively  a 
voltage  source. 


14.0  Boundary  Conditions 

In  practice,  one  knows  the  flux,  n+u,+  +n++u2  +  ttQvQz  This  is  not  conserved  (it  can  vary  if  the  thruster 

generates  nonsteady  flow)  but  it  is  constant  at  the  anode,  where  it  equals  the  gas  feed  rate.  If  ions  backstream  at 
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the  anode  (from  a  downstream  ionization  zone)  their  density,  velocity  and  temperature  are  obtained  from 
extrapolating  back  to  z=0.  Otherwise,  we  assume  that  the  incoming  gas  is  1%  ionized  and  0.1%  doubly  ionized. 


From  the  total  flux  and  the  ion  fluxes,  one  has  the  anode  neutral  flux.  The  neutral  density  at  z-0  is  obtained  by 
solving  simultaneously  for  density  and  velocity  to  match  that  flux.  Using  Eq.  [43]  for  the  neutral  velocity,  taking 
the  gradient  at  a  point  z}  by  differencing  the  density  at  points  zj+j  and  Zj.,,  and  assuming  that  the  gradient  at  z0=0  is 
the  same  as  that  at  Z/,  one  has  the  flux: 


Nf(Q)  =  no  (0)v0z  (0)  =  n0  (0)v0j 


«o(^2)-«o(°) 


2Az 

3*0  (OK 
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w  * 
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giving, 
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v0i  l  W  J  V  v0s 
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15.0  Initial  Conditions 

The  program  requires  initial  profiles  for  the  eight 
variables.  These  need  not  be  very  accurate,  but  they 
must  satisfy  the  boundary  conditions.  The  first  runs 
used  simple  constant,  exponential  or  Gaussian 
functions  with  a  few  variable  parameters.  The 
program  is  structured  so  that  the  results  of  one  run 
can  be  sued  as  the  starting  point  for  subsequent  runs. 
Although  the  program  evolves  the  variables  in  time, 
there  are  indications  that  the  initial  conditions 
matter,  that  the  discharge  has  more  than  one  mode  of 
operation  for  a  given  set  of  operating  parameters  and 
that  different  initial  conditions  can  lead  to  different 
steady-state  discharge  behaviors. 

16.  Solution  of  the  Model 

To  step  in  time,  the  program  starts  with  the  eight 
primary  variables  and  the  current,  calculates  E,  ne, 
Ve,  and  vfc  and  then  the  gradients  of  these,  as  well  as 
secondary  quantities,  such  as  all  the  collision  and 
ionization  rates.  To  make  the  calculation  stable,  the 


program  uses  either  forward  or  backward 
differencing,  depending  on  the  signs  of  the 
velocities.1221  Then  it  calculates  the  eight  F’s,  the 
diffusion  coefficients  and  finally  the  time  derivatives 
of  the  eight  primary  variables. 

The  program  also  integrates  E  (for  the  three  different 
currents)  to  obtain  the  discharge  voltage  and  then 
calculates  the  current,  J,  for  the  next  time  step.  To 
model  the  erosion  of  the  walls  by  sputtering,  it 
stores,  for  each  ion  species,  the  ion  wall  flux, 
0.6n+vL+cos(y/),  the  incident  velocity,  vL+,  and  the 
angle,  (y/+Q  that  the  incident  ions  make  with  the 
tangent  to  the  wall,  all  as  functions  of  z.  Finally,  it 
also  calculates  the  test  quantities  S )  and  5). 

We  have  begun  exploring  the  predictions  of  this 
model.  The  following  series  of  current  traces  show 
results  of  runs  in  which  the  voltage  was  varied  while 
the  other  parameters  (magnetic  field,  channel  width, 
gas  feed)  were  kept,  constant.  As  can  be  seen,  at 
lower  voltages,  the  discharge  exhibited  oscillatory 
behavior  while  higher  voltages  produced  a  steadier 
current. 
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Figure  3  200  V  discharge 
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Figure  4  250  V  discharge. 


Figure  7  400  V  discharge. 
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Figure  8  450  V  discharge. 


In  another  series  of  runs,  the  voltage  was  kept  constant  while  the  width  of  the  discharge  channel  was  varied.  The 
discharges  in  the  narrower  channels  approached  constant  current  while  those  in  the  wider  channels  exhibited 
oscillations. 


17.0  Conclusions 

A  detailed  comparison  with  experimental  data  has 
not  yet  been  made,  and  we  have  only  begun  to 
explore  the  parameter  space  accessible  through  this 
model,  but  the  results  obtained  so  far  are  in  general 
agreement  with  Hall  thruster  performance.  The 
success  of  the  derivation  described  here  and  the 
resulting  numerical  calculations  show  that  a  one¬ 
dimensional  Hall  thruster  model  can  include  much 
more  of  the  discharge  dynamics  than  had  been  done 
in  previous  1-D  models.  In  particular,  the  ion 
recombination  at  the  side  walls  and  the  effects  of  a 
varying  channel  width  and  a  non-radial  magnetic 
field  can  all  be  retained.  Although  this  adds 
complexity,  the  resulting  simulation  is  still  much 
simpler  and  therefore  runs  much  faster  than  a  full 
two-D  code. 


that  causes  sputtering  which  limits  the  operating 
lifetime.  It  stores  spatial  profiles  of  all  the  variables 
at  regular  times,  which  should  help  us  understand 
the  nonlinear  dynamics  of  these  devices. 
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